%% fig11b_phase_diagram.m
% 绘制 Fig.11(b): 在 x1-y1 初始平面下，不同初始条件对应的相图（以 x1 对 y1 为例）
clear; clc; close all;

%% 1. 设置系统参数（与论文一致）
params = struct(...
    'a1', 1.5, ...    % HR 神经元1 激活参数
    'b1', 3, ...      % HR 神经元1 非线性抑制系数
    'c1', 1, ...      % HR 神经元1 恢复变量
    'd1', 5, ...      % HR 神经元1 抑制作用
    'i1', 0.45, ...   % HR 神经元1 外部电流
    'a2', 0.77, ...   % FN 神经元 激活参数
    'b2', 0.3, ...    % FN 神经元 非线性抑制系数
    'c2', 0.8, ...    % FN 神经元 恢复变量
    'i2', 0, ...      % FN 神经元 外部电流
    'a3', 0.9, ...    % HR 神经元2 激活参数
    'b3', 3, ...      % HR 神经元2 非线性抑制系数
    'c3', 1, ...      % HR 神经元2 恢复变量
    'd3', 5, ...      % HR 神经元2 抑制作用
    'i3', 0.5, ...    % HR 神经元2 外部电流
    'm12', 0.785, ... % 耦合参数 FN -> HR1（或 HR1 ← FN）
    'm21', 0.52, ...  % 耦合参数 HR1 -> FN
    'm23', 0.785, ... % 耦合参数 FN -> HR2
    'm32', 0.08);     % 耦合参数 HR2 -> FN

% 固定其余初始条件（不在 x1-y1 平面变化的变量）
fixedVals = [3; 1; 1; 0.1; 0.3];  
% 解释：x2(0)=3, y2(0)=1, x3(0)=0, y3(0)=0.1, φ(0)=0.3

% 固定 FN 神经元参数 ε（论文中给定值）
eps_val = 4.87;

%% 2. 定义 x1(0)-y1(0) 初始平面上的不同初始条件
% 这里取三个代表性初始条件（可根据实际需要增减）
IC1 = [-2; 0];
IC2 = [2; 0];

ICs = {IC1, IC2};
colors = {'g', 'b'};  % 使用绿、蓝区分不同初始条件

%% 3. 积分参数设置
T_total = 5000;  % 总积分时间
tspan = [1000, T_total];
opts = odeset('RelTol', 1e-6, 'AbsTol', 1e-9);

%% 4. 对每个初始条件进行仿真，并绘制相图（取 x1 与 y1）
figure; hold on;
for k = 1:length(ICs)
    % 构造完整初始条件：X0 = [x1(0); y1(0); x2(0); y2(0); x3(0); y3(0); φ(0)]
    IC_xy = ICs{k};
    X0 = [IC_xy; fixedVals];
    
    % 使用 ode45 进行长时间积分
    [t, X] = ode45(@(t, X) hrfnhr(t, X, eps_val, params), tspan, X0, opts);
    
    % 绘制相图：这里以 x1 对 x2
    plot(X(:, 1), X(:, 3), 'Color', colors{k}, 'LineWidth', 1.5);
end

xlabel('$x_1$', 'Interpreter', 'latex', 'FontSize', 12);
ylabel('$x_2$', 'Interpreter', 'latex', 'FontSize', 12);
title('Phase Diagrams for Different Initial Conditions','FontSize', 14);
legend('IC1', 'IC2', 'Location', 'best');
grid on;
axis equal; axis tight;
